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Abstract. We review the current state of Cepheid modeling and discuss its dominant deficiency, namely the use of time 
dependent mixing length. Notwithstanding, Cepheid modeling has achieved some excellent successes, and we mention some 
of the most recent ones. Discrepancies between observations and modeling appear not so much in the gross properties of 
single mode Cepheids, but rather when more subtle nonlinear effects are important, such as in double mode or even triple 
mode pulsations. Finally we discuss what we consider the most important challenges for the next decade. These are, first, 
realistic multi dimensional modeling of convection in a pulsating environment, and second, the nonlinear modeling of the 
nonradial pulsations that have been observed, and, relatedly, of the Blazhko like phenomenon that has recently been observed 
in Cepheids. 
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INTRODUCTION 

When looking back to the pioneering days of Cepheid 
modeling {e.g., |H 01) we have to modestly admit that 
in terms of numerical modeling methods we have only 
made some progress, except that we now have at our 
disposal enormous computing power compared to that 
available in those days. 

There has however been quite a bit of progress in terms 
of the quality of the opacities that are now used (OPAL 
101, OP 0], and Andersen-Ferguson H]). In addition, 
some treatment of convection is necessary and is rou- 
tinely included in the calculations in the form of time- 
dependent mixing length. 

In parallel to direct numerical simulations, the ampli- 
tude equation formalism has been developed 0, H @> 
ilii (for a review see iflllp . The physical conditions 
that prevail in Cepheids and RR Lyrae, namely that the 
growth rates of the dominant modes are small compared 
to their frequencies, form the basis for the applicability 
of these techniques. We note in passing that, in contrast, 
amplitude equations do not apply in Pop. II Cepheids, RV 
Tau and Mira variables. In these stars the growth rates 
are comparable to the frequencies and therefore the am- 
plitudes can change substantially on the time scale of a 
period, hence the fre quen t occurrence of irregularity in 
the pulsations lfl2[ [Bj [l4|l . 

In many respects the application of amplitude equa- 
tions complements brute force modeling. Amplitude 
equations, or normal forms as they are known in non- 
linear dynamics lfl5ll . are very general and describe the 
underlying mathematical structure of the temporal be- 
havior, i.e., the pulsations in our case. In addition, when 
combined with numerical simulations iflrjll . they provide 
not only an excellent description of the modal selection 



problem, e.g., of where the region of double mode behav- 
ior occurs, or where hysteresis occurs ('either or' regime 
in the jargon of stellar pulsations), but they also yield a 
deeper understanding of these bifurcations in the pulsa- 
tional behavior lfl7ll . 



CRITIQUE OF THE HYDROCODES 

A number of hydrocodes have been developed over the 
last dozen years or so (e.g., the Italian code itla], the Vi- 
ennese code lfl9ll . the Florida-Budapest code [20] code. 
They are all essentially the same in that they include 
some similar form of time-dependent mixing length ap- 
proximation. The first code is based on the formulation 
of JH, and the second on that of JUIH, and the third 
a hybrid of the two formulations. Some of their fea- 
tures have been compared in ll24ll . More recent codes 
are the Australian I251 and the Polish ones l26ll . Some 
of these codes include a moving adaptive mesh that re- 
solves better the sharp temperature features, in particu- 
lar those associated with the hydrogen partial ionization 
front lfl9l l27ll . Some in addition improve on the equilib- 
rium heat diffusion J28ll . 

We briefly recall the ingredients of time dependent 
mixing length. The effects of turbulence and convec- 
tion in 1 -dimensional hydrodynamics description appear 
through the turbulent pressure and viscous stresses p t and 
p v , the convective F c and the energy coupling term c € . 
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The turbulent and convective quantities are assumed to 
be functions of the 'turbulent energy' e t that is assumed 
to satisfy 

^ + (A+*0£ = -A|-(^)-*' <3> 
at at pr L or 

The expressions for p t ,F t ,p v and ^ can be derived 
in analogy with gas kinetic theory, however without the 
same solid physical basis, or more simply just on dimen- 
sional grounds plus a few assumptions. All these terms 
contain (a) parameters of O(l) for whose values physics 
provides little guidance. There is also some ambiguity, 
e.g., possible physically acceptable choices for F c and the 
source term S, (that goes into the energy coupling term 

F c ~ a c e t ' 2 Y F c ~ a c e, Y [ ' 2 (4) 

S t ~a s e) l2 Y S,~a s e t Y l l 2 (5) 

where 

H„ dsi H n ds 

Y= or Y = -jr (6) 
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More a parameters appear if a Peclet correction for 
radiative losses and a flux limiter for F c < p c p T c smmc j 
are included. 

There are a number of well known problems with 
time dependent mixing length. To start, it is an empiri- 
cal rather than a consistent physical description. Second, 
the gradient approximations, i.e., F t « Ve t and F c °c Vs, 
are poor j29ll . In the same vein, 3D numerical simula- 
tions show that the upward and downward plumes are 
highly non local l30ll . whereas mixing length assumes 
a local approximation (diffusion of turbulent energy). 
Third, some 8 or more (a) parameters appear in time de- 
pendent mixing length and physics provides no guidance 
for their values. Therefore some of these parameters are 
calibrated with the help of a comparison of the results 
with observations, but the parameter space is large! The 
less important parameters are arbitrarily fixed. It is some- 
what disturbing that for best agreement with observations 
one needs a different set of a values for RR Lyrae and for 
Cepheids, and also for different metallicity. Finally, there 
is no guarantee that time dependent mixing length is flex- 
ible enough to describe convection sufficiently well, de- 
spite the many a parameters. 

But these problems notwithstanding dependent mixing 
length has achieved many successes. 

Before going on to discuss these successes we want to 
display in Fig. [TJthe behavior of the turbulent energy e t 
over a period in a pulsating Cepheid model with a period 
of 10.9 days, and M=6.1M , L=3377L , r eff =5207K, 
X=0.70, Z=0.02 l3lll . In the figure e t is shown as a func- 
tion of Lagrangean radius and the lightness of the grey 



reflects the strength of e t . The turbulent energy is largest 
in the region associated with the combined H and first He 
ionization fronts which is the convectively most unstable 
region, but it also appears in the next most important, 
viz. the He + -He ++ region. Fig. Q] clearly demonstrates 
how the turbulent energy tracks the source regions which 
move through the fluid during the pulsation. It also shows 
the importance of time dependence in the convective pul- 
sating envelope. The turbulent energy increases during 
the pulsational compression phase and that, in this rela- 
tively hot model, the two turbulent zones briefly merge. 

RECENT RESULTS 

The literature is very vast, so we will just mention some 
of the more recent work. 

One of the most striking properties of the Cepheids, 
both for fundamental (F) mode pulsations and for first 
overtone (Ol) mode pulsations is the progression of the 
Fourier decomposition coefficients, of the light curves 
as well as of the radial velocity curves {e.g., 1321 1. For 
F mode Cepheids this progression is of course related 
to the well-known Hertzsprung progression in the bump 
Cepheids. Full amplitude Cepheid model sequences do a 
good job at reproducing the Fourier properties of the F 
Cepheids JH and of the Ol Cepheids H. 

In parallel, the amplitude equation formalism that was 
developed for explaining the effect of internal resonances 
on the appearance of the light and radial velocity curves 
has given a clear demonstration that it is the 2:1 reso- 
nance between the self-excited F mode and the vibra- 
tionally stable, but resonant second overtone j3~5ll36ll that 
is responsible for the structure of the Fourier coefficients 
for periods around 10 days, rather than a Shockwave that 
reflects off the core. For the first overtone Cepheids (or 
s Cepheids) it is the 2:1 resonance of the stable fourth 
overtone with the self-excited first overtone that causes 
structure at period in the vicinity of 4 days. Modeling 
also reproduces well the observed ligh t and radial veloc- 
ity curves of individual stars, e.g., wfa , JUt]. 

A full amplitude model survey of F and Ol Cepheids 
1I39I1 has found that the phase lag between light and radial 
velocity curves is in good agreement with observations. 
This study has furthermore demonstrated that the phase 
lag can also be used observationally as a discriminant 
between F and Ol mode pulsation. 

The period ratio P\/Pq as a function of period Pq of the 
beat Cepheids is a sensitive function of the metallicity Z. 
This property has recently been taken advantage of to 
determine the metallicities in the LMC and SMC Hot 
with the help of Cepheid modeling fill l42ll . The same 
method has been applied to the 5 known beat Cepheids 
in M33. Interestingly, this yields a galactic metallicity 
gradient for M33 that is in good agreement with totally 





independent methods l43ll . 

A recent study shows that a comparison of the light 
curves of full amplitude bump Cepheid models with 
observed light curves can also be used as metallicity 
tracers and distance indicators iRill . 

Theoretical period-color-luminosity relations are in 
good agreement with observational ones l45ll . l46ll . A 
summary of the recent Italian work on this and other top- 
ics in Cepheid modeling appears in 14711 . 

Theory has been ahead of observations by predict- 
ing the existence of 'strange' Cepheids and RR Lyrae 
l48l l49ll . These are Cepheids or RR Lyrae that pulsate 
in a high (7th to 12th) overtone in which the pulsation is 
confined to the outer region, more specifically, above the 
hydrogen ionization front. Typically, the predicted peri- 
ods of these self excited modes are 4 to 5 times smaller 
than the fundamental period of the same object. The am- 
plitudes are predicted to be in the millimag range. On 
the observational side some evidence for the existence of 
strange Cepheids has been uncovered. However, it is hard 
to distinguish between intrinsic pulsation and ellipsoidal 
binary motion at the millimag level I50I [5TI0 . Further- 
more, in crowded field conditions, it is possible that the 
objects could be regular, albeit low amplitude Cepheids 
that appears above the P-L relation because of contami- 
nation by another star. Follow-up observations would be 
useful to ascertain that the identified objects are indeed 
strange Cepheids. 



In Fig. |2] we display the results of a recent search for 
ultra-low amplitude (ULA) variability in the combined 
MACHO and OGLE LMC data JH. For reference and 
to guide the eye, the well known Cepheid P-L strips are 
indicated as black dots. All the ULA (Fourier amplitude 
A < 0.01) objects that were found in the plotted period- 
Wesenheit magnitude (W = I — 1.55(V — /)) are shown 
as orange filled circles. It is possible that some of the ob- 
jects are ellipsoidal binaries that with the available data, 
are indistinguishable from ULA Cepheids. The cluster- 
ing of the ULA objects with the LMC Cepheids suggests 
though that the latter are indeed pulsating variables. The 
objects below the fundamental P-L relation are probably 
Pop. II Cepheids. The single object way above the P-L 
relation is a strange Cepheid candidate. When one con- 
centrates on the remaining objects, their position is seen 
to correlate strongly with the LMC Cepheid P-L relation. 
One expects any ULA Cepheids to lie at the edges of the 
instability strip, but Fig.|2]might suggest that they lie on 
a parallel P-L relation that is shifted slightly upwards. If 
that is indeed the case, it is an interesting challenge to 
explain the nature of these ULA Cepheids. 

The analyses of the OGLE data base e.g., SESEl), 
to which we refer the reader for the following discussion, 
have yielded a beautiful overview of the HR diagram. 
In addition to the F and Ol Cepheids, a whole zoo of 
other types of Cepheid pulsational behavior has now 
been observed, namely 02 single mode, double mode 
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FIGURE 2. Ultralow amplitude objects in the LMC, shown as open (red and orange) circles. The smaller (orange) circles possibly 
have proper motions. The classical Cepheids are superposed as dots for reference. 



F/Ol, double mode 01/02, triple mode, and Cepheids 
with Blazhko like modulations. 

For comparison we present the results of a recently 
computed HR diagram Isill in Fig. [3] The bounding 
curves on the right [in red, in the color figure of the elec- 
tronic version] and the [blue] ones on the left represent 
respectively the red and blue edges of the F, Ol, 02 and 
03 pulsational instability strips. We recall in passing that 
these boundaries are obtained from linear vibrational sta- 
bility analysis, and as such, they only tell us where at 
least one mode is linearly unstable, i.e., where pulsa- 
tions occur, but they generally do not tell us what type 
of pulsation actually results. For example, for F/Ol dou- 
ble mode behavior it is necessary, but not sufficient that 
both modes be linearly unstable (in the figure the inter- 
section of the F and the Ol instability strips). Nonlin- 
ear calculations are necessary to determine which pulsa- 
tional behavior can actually occur. (We say 'can occur' 
rather than 'occurs' because there may also be hystere- 
sis, i.e., the pulsational state depends on the evolutionary 
path.) 

Each of the symbols in Fig. [3] represents the result 
of one or more hydrodynamic integrations to steady full 



amplitude. The numbers on the left denote the mass (in 
Mq) which is constant along horizontal lines. The fun- 
damental periods (in days) are constant along the dotted 
lines. The open squares [black] denote the regime where 
only single mode F full amplitude pulsations (limit cy- 
cles) occur, the open upside-down triangles [red filled 
triangles] the regime of single mode Ol, and the open 
circles [blue filled circles] that of single mode 02 pulsa- 
tions. 

The solid [lightblue] squares denote the regime where 
either F or Ol pulsations occur, i.e., where there is hys- 
teresis. A star, say with mass 5M Q , evolving leftward en- 
ters the F instability strip (rightmost red line) and starts 
to pulsate in the F mode. It continues to do so to the left- 
most edge of the (lightblue) hysteretic region where it 
switches to Ol pulsations until it exits through the Ol 
blue edge on the left. On its return, moving now right- 
ward, it enters the Ol instability strip (leftmost blue line) 
and starts to pulsate in the Ol mode. It continues to do so 
to the rightmost edge of the (lightblue) hysteretic region 
where it switches to F pulsations until it exits through the 
F red edge on the right. 

Below, and as a continuation of the hysteretic region 
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FIGURE 3. Computed HR diagram indicating the domains of single mode pulsations, of double mode pulsations and of hysteresis 



we typically find F/Ol double mode pulsations denoted 
by crosses [yellow filled circles with a concentric black 
dot]. Care has to be exercised in the modeling compu- 
tations that the pulsations are truly double mode (i.e., 
with constant amplitudes and phases) rather than that the 
model may be switching from one mode to the other. In- 
deed, it has been shown that the amplitudes may level off 
relatively fast, and then vary extremely slowly, giving the 
erroneous impression of steady pulsations. To be certain 
that double mode is indeed obtained it is necessary to in- 
tegrate the model with several judiciously chosen initial 
kicks and to watch the evolution of the individual ampli- 
tudes in an amplitude-amplitude portrait lfl7ll . 

Interestingly, F/Ol double mode pulsations can also 
occur in a narrow region above the hysteresis regime. 
This double mode behavior is induced by the Pq = 2Pi 
resonance in the vicinity of the 10 d period. The locus of 
the center of this resonance is indicated as a dashed line 
in the figure. 

The width of the hysteresis regime [light blue lines] is 



clearly nonmonotone with mass. This is not a numerical 
problem, but we have again traced this behavior to a 
resonance, this time Pi = 2P4, as seen by the locus of 
the resonance center (solid line). 

Similar regimes of hysteresis and double mode occur 
on the right towards the lower masses, but now for Ol 
and 02. The hysteretic regime of Ol and 02 ('either 
01 or 02') is marked by solid triangles [green filled 
circles]. Below this occurs a regime of 01/02 double 
mode pulsations labeled with x's [orange filled circles]. 

We note in passing, and as a curiosity, that we have 
also found an extremely narrow strip of hysteresis with 
either F single mode pulsations or F/Ol double mode 
pulsations. 

It is likely that multimode and hysteresis also occur for 
the 2 Mo models, for which 03 is linearly unstable, but 
we have not checked this. 

Overall, the modeling is seen to reproduce the dif- 
ferent types of single mode and double mode behavior 
with approximately the right topography. However, some 



types of observed pulsational behavior, such as triple 
mode, have yet to be modeled. It is clear that double 
mode or triple mode behavior is much more sensitive to 
the details of the physical input than, say single mode, 
because it depends on more subtle effects, viz. the non- 
linear coupling between the modes. 

Before going on we also want to mention some well 
known, but serious discrepancies, for which stellar evolu- 
tion calculations are at fault, but on which Cepheid mod- 
eling unfortunately has to rely. 

First, low mass, low Z evolution loops do not pene- 
trate the instability strip where Cepheids are actually ob- 
served. It has been suggested that this could be a metal- 
licity selection effect |p5i [56ll . There remains disagree- 
ment about the treatment of convection and convective 
overshoots which leads to uncertainties in the Cepheid 
M-L-Z relation. Unfortunately, this in turn causes uncer- 
tainties in Cepheid modeling which depends on an ML 
relation. However, we note that the modeling of LMC 
and SMC Cepheids lf57ll based on the tracks of Girardi 
et al. ll58ll is in good agreement with the resonance con- 
straints imposed by observations. Finally, the occurrence 
of a bump mass discrepancy was largely removed with 
the OPAL opacities |5')|. some authors claim that there 
remains a small discrepancy |44]Lwhich is at variance 
with the above mentioned work 15711 . 



BEYOND TIME-DEPENDENT MIXING 
LENGTH 

Since dependent mixing length has such a long list of 
known problems, efforts have been made to go beyond 
the simple model described above, while keeping its ID 
feature. This unfortunately leads to a proliferation of 
equations |29[ l60l l6ll l62ll with a concomitant numerical 
cost and perhaps little guarantee of substantial improve- 
ment. 

In parallel, plume models have been proposed l63ll . 
Along similar lines Stoekl has developed a 2 column 
model l64ll . However, none of these approaches have 
been used in Cepheid models, nor are they expected to 
provide a parameter free description of convection which 
after all is 3 dimensional phenomenon. 

TWO GRAND NUMERICAL 
CHALLENGES 

In our opinion there remain two grand challenges for the 
modeling of Cepheids, and for that matter of RR Lyrae 
which are quite similar in many ways. These are the mod- 
eling of convection in the highly structured envelopes of 



Cepheids, on the one hand, and of nonradial pulsations, 
on the other. 

1. Multidimensional simulations of 
convection in Cepheids and RR Lyrae 

There is a great deal of literature on the numerical 
difficulties that arise in the modeling of astrophysical 
convection. The literally astronomically large Rayleigh 
numbers imply that the turbulence is well developed. 
This in turn implies that many degrees of freedom are 
involved which sets extreme mesh requirements. How- 
ever, the hope of all large eddy simulations (LES) is that 
once we have attained a sufficient spatial resolution, then 
the dominant, large scale features of convective transport 
are well reproduced. Another problem arises from the ex- 
tremely small Prandtl number. (The Pr number is the ra- 
tio of gas viscosity to heat conduction.) In the numerical 
modeling the viscosity is however not set by the physical 
extremely small viscosity, but by that of the numerical 
scheme. Effective Pr numbers of order 0.01 or less are 
hard to achieve. 

There has been a great deal of excellent multi- 
dimensional numerical hydrodynamics work in other as- 
trophysical contexts, such as solar models, stellar core 
burning, dynamos 11651 16611 . Additional difficulties arise 
that are specific to Cepheids and RR Lyrae modeling. 
First, it is out of question to compute convection over 
the whole sphere, nor is it probably necessary. If con- 
vection occurs over a vertical range AR C it seems reason- 
able that we limit ourselves to a horizontal sector whose 
width is a few times AR C at the vertical center R c of the 
convective region, i.e., an angular sector AO = AR C /R C , 
typically a degree or so, and impose periodic horizontal 
boundary conditions. But even this leaves extreme mesh 
requirements. The resolution of the rapidly varying tem- 
perature structure in the partial hydrogen ionization re- 
gion requires tiny AR. However, the hot inner region that 
needs to be included in a realistic stellar model, even 
though no convection is taking place there requires large 
AR because the sound speed is large, so that we are not 
killed by the Courant condition. The mesh aspect ratios 
can thus vary from 0.01 to 100. 

Second, the surface boundary condition is delicate be- 
cause the envelope is surrounded by an optically thin, ap- 
proximately isothermal atmosphere with an exponential 
decrease of density. 

Third, the pulsating environment in which the stellar 
radius can change by more than 10%, combined with the 
exponentially decreasing atmosphere, make an Eulerian 
code unsuitable. A moving mesh is needed to follow 
the average pulsational motion, and adaptive features 
are required to track the sharp temperature gradients (H 



ionization front). 

It should be clear that the modeling of convection in a 
pulsating star is a challenging problem that is at the limit 
of current modeling possibilities. However an intermedi- 
ate goal could be to make a small but sufficient number 
of 3D simulations in realistic Cepheid or RR Lyrae mod- 
els, and then use the results to calibrate the mixing length 
recipe parameters, or if necessary to improve or even re- 
place mixing length theory. 

2. The Modeling of Nonradial Pulsations in 
Cepheids 

The analysis of the observations of the OGLE 
Cepheids strongly suggest the excitation of nonradial 
modes HHH. 

On the theoretical side, the linear stability analysis of 
nonradial modes in Cepheids and RR Lyrae is notori- 
ously difficult Ic37ll . Ic38ll . The problem arises from the 
fact that the nonradial modes have mixed p and g mode 
nature with unresolvably rapid spatial oscillations in in- 
terior. This means that most nonradial modes are very 
strongly damped. However, some surface-trapped modes 
were found to be unstable or only mildly damped, and 
they can therefore compete dynamically with the low or- 
der radial modes ll69l F70h . A recent revisit of the prob- 
lem with realistic stellar models lf7lll however finds no 
nonradial instability among the low I modes. (Only low 
I are detectable with full disk observations.) This prob- 
lem therefore remains a puzzle even at the linear level. 
Ultimately, of course, nonlinear 3D hydrodynamic simu- 
lations will be necessary to model the observed nonradial 
pulsations, which is again a very tough numerical prob- 
lem, for many of the same reasons as for convection 

The analyses of the LMC variables f52ll53ll also show 
the existence of Blazhko like modulations of the pulsa- 
tions which bear some similarities with the same effect 
in RR Lyrae. Just as in the latter the physical origins re- 
main obscure. Similar modulations also seem to occur in 
in V473 Lyr and in Polaris ||72|l . 

Among the viable models for the Blazhko effect are 
the likely involvement of one or more nonradial modes 
lf73ll . A specific mechanism involving resonant coupling 
has been examined by |p74tl . 

One may also wonder whether the nonradial modes 
can be destabilized by convection? Could this be the 
origin of the observed Blazhko modulations f53ll ? We 
are somewhat skeptical that a simple dependent mix- 
ing length description is sufficiently good to answer that 
question. Most likely a full 3D simulation, as described 
in the previous section, will be necessary to treat the con- 
vection - pulsation interaction. There is also the possibil- 
ity the Blazhko modulations arise from the radial pulsa- 



tion - convection interaction. In any case an explanation 
of the physical mechanism of the Blazhko cycle remains 
a theoretical challenge. 

Despite these serious difficulties, we think that these 
two grand modeling challenges will be attacked and pos- 
sibly will be met in the next decade. 
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